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The one dimensional Coulomb lattice fluid in a capacitor configuration is studied. The model is 
formally exactly soluble via a transfer operator method within a field theoretic representation of 
■ the model. The only interactions present in the model are the one dimensional Coulomb interaction 

between cations and anions and the steric interaction imposed by restricting the maximal occupancy 
at any lattice site to one particle. Despite the simplicity of the model, a wide range of intriguing 
£NJ ■ physical phenomena arise, some of which are strongly reminiscent of those seen in experiments and 

numerical simulations of three dimensional ionic liquid based capacitors. Notably we find regimes 
where over-screening and density oscillations are seen near the capacitor plates. The capacitance is 
also shown to exhibit strong oscillations as a function of applied voltage. It is also shown that the 
corresponding mean field theory misses most of these effects. The analytical results are confirmed 
by extensive numerical simulations. 
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I. INTRODUCTION 
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The physics of interacting Coulomb systems has been extensively studied for decades [IH1]- Most of this study 
has been concentrated on the statistical mechanics of electrolytes such as salt solutions. While there still exists a 
number of outstanding questions concerning the bulk thermodynamics of electrolytes, there has been just as intense 
an effort to understand how the presence of electrolytes modify the interaction between charged surfaces, colloids and 
macromolecules. The underlying statistical mechanics problem being to date unsolved, the study of these systems 
relies on the use of approximations. The simplest models of electrolytes are called primitive models, where the ions 
are taken to be point-like charges and the only interactions are electrostatic. 

For these simple electrolyte models with charged interfaces, two regimes are now reasonably well understood and 
their results have been validated numerically. The first is the so called weak coupling regime where the mean field or 
Poisson-Boltzmann (PB) theory [4j correctly describes the system, ignoring fluctuations about the mean field. These 
can be taken into account to refine the accuracy of the results. The PB theory is valid when the distance over which 
the ions interact with each other (the Bjerrum length) is smaller than the distance at which they interact with the 
charged surface (the Gouy-Chapman length). From another mathematical point of view the field theory describing 
interacting charged particles can be analyzed in the saddle point approximation pj, and in the regime of its validity 
the fluctuations renormalize the mean field results improving its accuracy. The other limiting regime is the so called 
strong coupling regime where the problem can be treated within an effectively virial-like expansion based upon the 
^ interaction of single particles with the charged surface @. This regime occurs when the Bjerrum length is much 
large than the Gouy-Chapman length, which becomes the case as the valency of the ions increases. The underlying 
instability in a point-like charge model due to the pairing of anions and cations means that the strong coupling 
approximation has only been applied to the counter-ion only case in the presence of charged surfaces. The region 
between these two limiting regimes is less well understood and a theory extrapolating between the two is still lacking. 

Results of calculations in these two limiting regimes do not have to be made finite by invoking a finite size of the 
ions. This is because (i) point charges are smeared out by the PB approximation and (ii) as stated above, the strong 
coupling approximation pertaining to the counter-ion only case, is effectively a one particle theory. However, it is clear 
that in certain circumstances steric effects due to the finite size of ions will come into play, not only especially near 
highly charged bounding surfaces but also in the bulk for systems composed of large charged molecules such as room 
temperature ionic liquids (RTILs). The presence of a new length scale given by the molecular size now complicates 
the above discussion of strong and weak coupling regimes, rendering the study of RTILs much more complex than 
simple electrolyte models. 

RTILs are fluids with large, asymmetric ions where steric effects cannot be neglected. These liquids have a number 
of practical technological applications @; for instance, they are used as electrolytes in fuel and solar cells as well 
as super-capacitors. The study of bulk properties of ionic liquids was stimulated by attempts to understand critical 
behavior of electrolytes. Here the effect of finite size ions needs to be taken into account in order to regularize the 
model. The restricted primitive model incorporates hard core effects in a continuum context and can be studied via 
approximate methods from liquid theory, for instance integral equations as well as field theoretic approaches (8l-[l0|. 
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Another way to take into account the finite size of ions is to restrict ions to points on a regular lattice (for a historical 
sketch of lattice Coulomb fluids (LCF) see [ill, [HI) and forbid that there be more than one particle at any given site 
[UHlfl]. This is clearly a rather idealized model since it assumes an underlying lattice structure which will influence 
the thermodynamics of the system and also since it treats the effective molecular sizes of cations and anions as the 
same. 

The study of LCFs in the capacitor configuration is a more recent subject of study. It is physically clear (T3] that 
the size of the ionic species will restrict the maximal local charge density that can be locally attained in a RTIL, for 
instance the counter-ion charge density near a highly charged surface will have a maximal value at closest packing. In 
certain regimes therefore, we should expect a fundamentally different behavior of ionic liquids at charged interfaces 
compared with aqueous electrolytes. Not withstanding this proviso, as pec ts of the behavior of RTIL capacitors can be 
captured by a mean-field treatment of the lattice Coulomb fluid model [l7l [l8| . Qualitatively the mean-field treatment 
of the LCF predicts that steric constraints cause the capacitance of ionic liquids to decay at large applied voltages 
(owing to the saturation of counter-ions at the capacitor plates). At the point of zero charge (PZC) the capacitance 
can exhibit a maximum and then decay monotonically, however in certain cases it can also be a local minimum, 
increasing to a maximum value at finite charge and then decaying (and is thus a non-monotonic function of applied 
voltage). Which of these two capacitance behaviors is predicted by the mean field theory depends on the lattice 
packing fraction [l9j]. For dilute electrolytes the capacitance at PZC is always a minimum (l7) . steric effects which 
reduce the capacitance only becoming important at higher voltages. The basic LCF model also uses a rather simple 
description of the capacitor plates, neglecting dispersion interactions that will come into play for ions near metallic 
electrodes and will influence the behavior of ions even at zero applied voltage [2(J HH ■ Taking into account these 
dispersion effects, that essentially lead to binding of discrete ions to their image charges in the metallic plate, leads 
to significantly different predictions for ionic distributions at the capacitor plates and consequently to significantly 
enhanced capacitances. The effects of these interactions improve the quantitative agreement between theory and 
experiment. 

At an experimental level, X-ray reflectivity [22j], SFA [23[ and AFM [24] studies at charged interfaces show an 
alternating charge distribution starting with an over- screening cationic layer, at the negatively charged substrate, 
which decays roughly exponentially into the bulk liquid with a periodicity comparable with the size of ionic species 
[25j . Similar results are found in numerical simulations [lj| . 

The charge layering seen in experiments and simulations is expected to be a generic feature of RTILs at sufficiently 
strongly charged interfaces resulting from an interplay between steric effects and strong electrostatic coupling. While 
the phenomenology of the differential capacity is described by the mean-field solution of the LCF model [ijj, the 
over-screening and alternate charge layering at a capacitor plate are not. 

In order to explain overcharging and charge oscillations, but staying within the mean field approximation, Bazant et 
al. [26[ proposed a phenomenological theory based on a Landau- Ginzburg-like functional containing the standard LCF 
free energy [27| but with an additional higher order potential-gradient term, similar to what is found in Cahn-Hilliard 
models. The origin of this higher order gradient term can be justified [28[ via a decomposition of the Coulomb 
interaction into a long-distance mean- field-like component and a non-mean-field strong coupling component [29j. 
This alternative interpretation of the origin of this phenomenological term as a result of strong coupling or correlation 
effects leads to the natural question as to whether it is the model or the approximation that needs to be modified to 
understand the physics of these systems. The latter possibility suggests that a more detailed non-mean-field analysis 
of the original Kornyshev LCF model is appropriate. Therefore, in this paper, we will study the exact solution of the 
one dimensional (ID) Coulomb fluid or gas where all the physics of the basic LCF model can be analyzed. While the 
physics of this system will be different to the 3D system, it is of interest because we can assess the validity of mean 
field theory in this case by comparing it with an exact solution. In addition the ID LCF corresponds to sheets of 
charge (at discrete lattice sites), and as experimental results and simulations suggest that a layering phenomena is 
present, one might hope that this model does contain and possibly explain some of the phenomena of the original 3D 
system. Indeed we shall see that the phenomena of over-screening, charge oscillations and some of the more exotic 
behavior of the capacitance of these systems are predicted by our exact solution. The conclusion of the paper is that 
it is quite possible that some of the physics of RTIL capacitors can be explained by the LCF model if the model is 
properly treated at the mathematical level. 

We thus propose an exact analysis, albeit in ID, that demonstrates the full physical phenomenology of the LCF 
model, and in particular show that charge density oscillations, over-screening, and oscillations in the capacitance 
emerge naturally from this model without the introduction of any new physical interactions. A short version of this 
work can be found in 30], the version here fully details the calculations as well as presents a new mathematical 
analysis of the exact solution, in particular to explain the surface charging through a fixed applied voltage for large 
capacitor plates. The physical properties of the systems are discussed in greater detail as well as the relationship 
between the exact result and the corresponding mean-field theory. Finally, given the rather rich and exotic behavior 
found in our analytical approach we have confirmed our analytic results by comparing them with extensive Monte 
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FIG. 1: Lattice Coulomb fluid model. Dark balls stand for cations and light balls for anions. Boundary charges are at sites — 1 
and M. 

Carlo simulations of the ID LCF. 



II. ONE DIMENSIONAL LATTICE COULOMB FLUID MODEL 

Here we describe the one dimensional LCF and show how it can be formally solved exactly. The first exact solutions 
of the one dimensional Coulomb fluid or gas are due to Lenard and Edwards [3l| where the bulk properties of a point- 
like model of ions was studied. The problem was then analyzed in a more mathematical context [32[ , and the peculiar 
dependence of the bulk free energy on the electric field applied to such a system was elucidated. Later the effective 
interactions between charged surfaces in this model was studied [HI, [34j in the presence of counter and co-ions and 
counter-ions only. These studies were also used to benchmark the range of validity of the various approximation 
schemes (strong and weak coupling) discussed in the Introduction. Here we present a formally exact solution to the 
LCF in one dimension based on a functional transfer matrix formalism which can be used to extract, numerically, all 
thermodynamical observables to any required numerical precision. 



A. Without external charges 



We consider a one dimensional system where sheets of density <7* (cations) or — er* (anions) can sit at points on a 
lattice of size M and lattice spacing a. If we assume that the system has an effective surface area A, the effective one 
dimensional Coulomb Hamiltonian can be written as 



H = -'- 



Aal M ~ X 



4e 



i,j=0 



■j\SiSj 



(1) 



where e is the dielectric constant within the bulk of the ionic liquid. 

In the above, Si is a classical spin taking the value Si = 1 if there is a cation at lattice site i, Si = — 1 if there is an 
anion and Si = if it is empty. The basic ID model is shown in Fig. [TJ We will impose overall electroneutrality on 
the system which means that Si = using the discrete Fourier representation 



5 £ ; Si,o 



1 
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dtp exp tip 




We will work in the grand canonical ensemble where the grand partition function can be written as 

Xfio 1 is* 



S = Tr 



s, 




\i - j\SiSj + iip > '.Si | r 2 - 



M-l 

E 

i=0 



2^ 



(2) 



(3) 



where fi is the fugacity of both types of ions (as the liquid is symmetric between anions and cations) . We have also 
introduced the dimensionless parameter 
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2e ' 



(4) 



which is the ratio of the electrostatic energy between two ions at neighboring sites to the typical thermal energy 
ksT = 1//3. Note that the overall charge of an ion in this notation is q = Ao* and in this notation 



_ q 2 /3a 
7 ~ ~2l7' 



(5) 
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Above the trace is over the spin values Si at each site and the integral over the field ip enforces the electroneutrality 
constraint. We notice that the integral over -0 may be translated, it just has to run over an interval of length 2n. 
Carrying out a Hubbard- Stratonovitch transformation gives 



CX P ( | E N - JlSiSj \=nJ [d4] exp i-^j 



ij=0 





d(f>(x) 


27/ 


dx 



dx + iS^Sj(j>{j) 



(6) 



where N is a normalization constant which can be determined at the end of the computation imposing that the grand 
partition function without ions (/i = 0) should be equal to one. The field 4> is related to the electrostatic potential of 
the problem. The energy of a point charge % in an electrostatic potential V{x) is simply given by £es = <}V(x) which 
for this system gives a factor in the exponential Botlzmann weight —(3£es = V(i)qi, comparing with Eq. © 

now allows us to identify up to an overall constant C, the electrostatic potential as 



Notably the potential drop across the system is given by 

AV = V(-l) - V(M) 



(7) 



(8) 



where we include the sites —1 and M which represent the capacitor plates. 

With electroneutrality, the integrand is invariant by a translation of the field <p(x) — > <j>(x) + 4>o- To keep a finite 
expression, we must set this constant enforcing, say, 0(0) = 0. Inserting this expression in the grand partition function, 
we can easily perform the sum over the Sj which are no longer coupled: 



(9) 



For this path integral, where the field <j> propagates freely between the integers, the propagator is known exactly 



[#]exp 



d<f>(x) 



dx 



/27T7 



(4>{j + 1) - d>{j)f 
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This now gives 



f TT d ^3) 

u 2^ 



(0(i + 1) - 
27 



^log(l + 2^cos(0(j) +tl>)) 



Finally we can rescale the variables <fi(i) + ip — > <j>{i) in the above integral to give 



(10) 



(11) 



d<t>{j) 



2tt I 

3 



exp 



E 



w + 1) - 

2 7 



^log(l + 2/icos(<K?))) 



(12) 
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All integrals here are for 4>(i) for i > 1 over the interval (—00, 00) but the integral over 0(0) is over (— 7T, 7t). 



B. With charges on the boundaries or constant applied voltage 



Charges can be put directly on the boundaries of our system by inserting a charge qQ on the site —1 and — qQ 
on the site M. To obtain the partition function with these charges, we just have to add iQ(<f>(— 1) — <j>(M)) in the 
exponential in (|12[) and extend the integration to over <fi(—l) and <fi(M), where we again restrain the integral over the 
field at the first site (<j>(— 1)) which is now at i = —1 over (— 7r, n) . This gives the grand partition function at fixed 
external charges 

*Q = J II tS= cx p (" E W + ^ 0O ' ))2 + E ^ + 2 ^ + - W) j • (13) 
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The grand potential is defined as 

r» Q = -/3- 1 log(S Q )=r 1 WQ , (14) 

with ujq being the dimensionless grand potential. 

In this capacitor set up we will study the differential capacitance 

c = "§v- < 15 > 

and for simplicity we will work with dimensionless quantities: the reduced potential 

v = f3qV = (16) 

and the reduced capacitance 

dQ 



dAv' 



(17) 



We can also impose a potential drop and look at its effect on the charges on the boundaries. To obtain the 
corresponding partition function, we have to restrict the integration in Eq. (I12[) only to be over fields satisfying 

Ay= WM) -«-!)) 

pq 

or in terms of the dimensionless quantities 

cf>(-l)-<j>(M)=iAv. (19) 
We now note that this constraint can be imposed via a delta function written as a Fourier integral 

<5(<H-1) - <P{M) - iAv) = -!- / dQ exp(iQO(-l) - 0(M) - iAv). (20) 

27T J 

The partition function and the grand potential at fixed potential difference can then be written, up to an unimportant 
factor of 2ir, as 



E A „ = exp(-u Av ) = / dQ exp(QAv -uiq). 



(21) 



The interpretation of this equation is clear, the system is now free to accumulate charge at the plates in order to 
hold the potential constant. This result is in accordance with standard thermodynamic relations and is exploited in 
numerical simulations of systems with fixed voltage drops across them [35j. 

There is a clear difference between the two ensembles. In the constant charge ensemble the potential drop Av is a 
fluctuating quantity and from Eq. (1191) its average value is given as 

(Av) Q = H(0(-1) - J>(M))) = (22) 

However in the constant potential ensemble the surface charge Q is a fluctuating quantity whose average value is 
given by 

(Q)a v = -^- (23) 
oAv 

For each ensemble we have a corresponding differential capacitance. At fixed Q we have 

d{Av)Q 



and at fixed Av we have 



d(Q)Av ,„-, 

cav = o~X • (25) 

oAv 
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In this latter constant potential ensemble we have a simple fluctuation dissipation relation relating the capacitance 
at constant voltage to the fluctuations of the charge on the capacitor plates: 

CAv = (Q 2 )av - (Q)av > 0- (26) 

We should note here that one of the first explanations of double layer capacitance is due to Helmholtz [38| who 
assumed that the surface charge could be neutralized exactly by a layer of ions at a distance r from the surface, where 
r is the effective ionic radius of the ions (assumed to be the same for both cations and anions). Applying this to our 
model (note the complete neutralization will not always be possible) we find a dimensionless Helmholtz capacitance 

c « = ^ ( 27 ) 

where in this non-fluctuating case cq and cav are the same. This Helmholtz capacitance sets an intrinsic scale of 
capacitance in the model studied here. 



III. MEAN FIELD APPROXIMATION 



In the mean field approximation, we compute the partition function by the saddle-point approximation, minimizing 
the argument of the exponential function. In Eq. (|12[) , this gives the following equation for the potential (we write 
vj = v(j)): 

2 Vj - v j+1 - t>j_i 2^sinhfa) _ Q 
7 1 + 2/icosh(wj) 

The reduced grand potential (measured in units of ksT) is then given as 



u> = - ln[H] = {V ' + \ V])2 - E 1 °g( 1 + 2^cosh( Wj )). (29) 

j j 

When the boundaries carry no charge or the applied potential Av = 0, the solution to this equation is obvious: vj = 0. 
Then the pressure is easily obtained 

p = /3P = -|^=log(l + 2 M ), (30) 

where the above is a discrete derivative and the formula is valid in both the constant Q and constant Av ensembles. 
In this expression, we note that the charge carried by the particles does not enter in the pressure which is, of course, 
due to the mean-field approximation in the bulk. 

In the case of a constant applied potential the boundary condition is simply given by v—i — vm — Av. In the 
presence of charges on the boundary the fact that the system is globally electroneutral means that the electric field 
is zero outside the system giving vq — u_i = vm — «m-i = lQ- 

Within the mean field approximation the local the charge density is given by 

2£smhfa) 

1 + 2/xcoshfe)- [6) 



There is a slight difference between this mean field approximation and that found in [17|, [19|, [27| where the continuous 
limit a — > with the maximal density q/a = po kept constant is taken. In this continuum approximation, in dimension- 
full quantities the mean field equation becomes 

d 2 V(x) 2 f i S mh((3qV(x)) 
6 dx 2 PQ l + 2ficosh((3qV(x))' [ ' 

This limit is legitimate if the potential varies slowly compared to the length a, but is not correct in the general case 
[27j . In order to compare the mean field theory with the full underlying model we will numerically solve the finite 
difference equations to make as thorough a comparison possible between the mean field theory and our exact solution. 



IV. EXACT COMPUTATION 



In what follows in order to lighten the notation we will denote yt = <j>{i). We define the operator 

i / (y-i 

— — cxp — 

/2Fy V 2 7 



p(y,y') ^=ex,.f -V^-l (■»■>» 



and also its operator square root 

p v W) = - 

V 7 
which obeys 

| dyV /2 (y, 2/V /2 (j/, y") - p(i/, y"). (35) 
With these definitions the grand partition function in the constant Q ensemble can be written as 

„ M 

5q = / dy- 1 (l + 2ij,cos(yo))Y[dy j p(y j -i,y j )(l + 2ij,cos(y j )),exp(iQ(y-i - y M )) (36) 
J 1=0 

where the integral over is over (— ir,ir) and the other integrals are over (—00,00). Mathematically it is now 
convenient to introduce the symmetric operator 



K{y,y') = J P 1 / 2 (y,z)(l + 2^cos(z))p 1 / 2 (z,y')dz 7 (37) 

using which we may write 

E Q = — [ dxaq>(iQx) \p 1 / 2 K M p 1 / 2 } exp(-iQar), (38) 
where we have used operator notation of the form 

/oo 
0(y,y')f(y')dy', (39) 
-00 

for the operators O which are K and p 1 / 2 in the above. Now, if we define the function 

*q{x) = -^=P 1/2 expHQx), (40) 

we can write 

~ Q = ^ dx V Q (x) [K M V Q ] (x) ee (* Q \K M \* Q ). (41) 

The computation of the local dimensionless charge at the site i is carried out as follows. In the carrying out the 
trace of the spins once decoupled we have 

Siexp(i(j>iSi) =2ishM>;), (42) 

S<=±1 



and we thus find that 



where 



_ (9 Q \K*LK"-*-i\9 Q ) 



L(y, y') = 2iii J dz p 1 / 2 {y 1 z) sin(s)p 1/2 (s, y'). (44) 
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For practical computations it is best to analyze the action of the various operators K and L on functions decomposed 
in term of Fourier modes. We define a Q-dcpendent class of functions since the vector space of functions having the 
Fourier representation 

f(x)= /fcexp(ifcx)- (45) 

kez-Q 

We see that the classes where Q have the same non-integer part are the same. The function is of Q class as 



p 1 / 2 exp(— iQx) = cxp ^— cx p(~ iQx), 



and thus 



*Q,fc = 



exp - 



">k,Q- 



In this representation we then have 



Kf{x) = ^ Kf k exp(ikx), 



(46) 



(47) 



(48) 



with 



Kf k = cxp(- 7 fc 2 /4) exp(- 7 fc 2 /4)/ fc + fi 



exp( )f k -i + cxp( )/ fc+ i 



Similarly in the same representation we have for the charge density operator L 

r f , ^ k \( , 7(fc-i) 2 w - , 7(fc + i) 2 w - 

i/fc = A*cxp( — ) cxp( )/ fe _i - cxp( )/ fc+ i 



In terms of matrix elements we can then write 



K k ,i = cxp( ) (5k,i + fJ-[Sk,i-i +$k,i+i\) , 



L k j = /xexp(- 
The grand partition function now reads 



4 

-f[k 2 + I 2 ] 



){5k,i+\ - Sk,i-i)- 



s Q = exp 



(c Q \K M \c Q ), 



(49) 
(50) 

(51) 
(52) 

(53) 



where the surface charge vector \cq) has components CQ.k — ^Q,k- The local charge density at site i is then given by 

(c Q \K*LK M -*-i\c Q ) 



{c Q \K M \c Q ) 



(54) 



The above expressions are formally exact and are given in terms of infinite dimensional matrices. In practice we 
now need to numerically evaluate the expressions in order to obtain the exact thermodynamics of the model. In 
practice we see that modes with large n are exponentially suppressed and thus we can truncate the matrices involved 
at sufficiently large n. 



A. Other formulation and limit of highly charged boundaries 



From another point of view the representation above presents a way of examining the problem in the strong coupling 
limit as a discrete path integral. The matrix products involved can be expressed over a discrete path integral of a 
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discrete random walker X n which at each step changes by AX n = ±1 or AX n — 0. This corresponds respectively to 
the presence of cations, anions or holes. The grand potential at fixed charge Q can be written as 

S Q = ]T exp(-^^(X„-Q) 2 + ln(/i)^(X„-X„ +1 ) 2 ) , (55) 

X-i=0,X M =0 \ n n / 

where the initial and final points of the process X n are fixed at zero. This representation is essentially a ID Coulomb 
lattice gas version of the Villain representation for the 2D Coulomb gases [36|, H3] and is also related to the electric 
field representation of the Coulomb gas used by Aizenman and Frohlich [32| . In the limit of large 7 it is clear that the 
process X n will be dominated by paths which start at zero and then stay as close as possible to the value Q. If we 
define by 9 the maximal distance that Q is from an integer we can write Q = uq + 9 where by definition 9 G [—1/2, 1/2] 
and riQ is integer. For large 7 the process X n will increase from X_i = to X nQ -\ = nq and similarly descend from 
this value when it is at distance uq from the end of the path m — N. These initial and final parts of these paths have 
a contribution to the action S of Sq (where Sq = exp(S')) 

tiq riQ 

Stait = -7^(Q-«) 2 + 21n(/i)n Q = -7^(71 + \9\ f + 21n(/x)n Q . (56) 

n=0 n=0 

This first term is independent of M and is not extensive in the system size, it can thus be identified with a surface 
term in the free energy. Now, in the bulk there are two possibilities: the first is where 7 3> ln(/^) — 1 and in this case 
it is not energetically favorable to depart from the optimal value uq, the effective bulk action thus being 

S h ^ = ~{M-2n Q )0 2 . (57) 

We thus see that in this limit the bulk free energy depends strongly on the surface charge. In the limit of large 7 
but where ln(/x) 3> ^9 2 , the system reduces its action by permitting charge oscillations about the optimal value tiq, 
giving a bulk action of 

Sbuik « -~(M - 2n Q ) (9 2 + (1 - \9\) 2 ] + (M - 2n Q )ln(/x). (58) 

This unusual dependence of the bulk thermodynamics on surface terms can be seen in an alternative manner. In the 
limit of large M the thermodynamics in the constant Q ensemble is dominated by the eigenvector |Ao(<5)) of maximal 
eigenvalue A (Q) of the operator K (thus K\X (Q)) — X {Q)\Xo(Q)))- The partition function in the constant Q 
ensemble obeys in this limit 

ln(S Q ) = -T±. + Mln(A (Q)) + ln((A (Q)|c Q ) 2 ). (59) 

The subtle point in Eq. (I5U1) is that the largest eigenvalue and eigenvector actually depend on Q. This phenomenon 
occurs in the ordinary continuum one dimensional Coulomb fluid without hard core interactions, the thermodynamic 
limit thus depends on the boundary term. Physically this is related to the long range nature of the one dimensional 
Coulomb interaction and in terms of field theory the ground states or 9- vacua depend on the boundary terms [32| . 
The system is thus sensitive to the ability of the counter-ions to completely screen the surface charges and from our 
discussion above the eigenvectors and eigenvalues of K depend on the non- integer part of Q. If one is in the strong 
coupling regime one may truncate the transfer matrix K by keeping just the site uq and its neighbors. This gives the 
truncated matrix 

/exp(-i(l + 2|0|)) /iexp (-2(1 +2|0|)) \ 

^ t =cxp(-^ 2 ) M exp(-2(1 + 2|0|)) 1 Mexp(-2(1-2|0|)) . (60) 

2 V ^exp(-2(l-2|0|)) exp(-2(l-2|0|))y 

In the limit of finite fi but large 7, and when 9^0 (more precisely when ln(/i) — 1 <C j0), we see that the matrix K 
can be approximated as 












1 


• 










X t «exp [-j-9 z \ 1 , (61) 
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and, in agreement with our argument above about dominant paths, we see that in this regime we have 

In(A (Q)) 



2 



(62) 



This means that at relatively small values of \i integer values of the surface charge are thermodynamically preferred 
and thus in the constant potential ensemble we should see plateaus of roughly constant integer charge in the average 
surface charge as a function of applied voltage. However, still in the strong coupling limit of large 7 but close to full 
filling where ^texp(— 7/2) 3> 1, we have 



K, -,,,xj»(-^ 2 ) | M ex P {-l(l + 2\6\)) 







and for non-zero 9 we thus obtain 



ln(A (Q)) 



/1 exp l 
/xexp 1 

\0\- : 



-i(l+2|*|)) 

/x exp (-5(1 -2^))), 

-i(l-2|0|)) 



ln(pi), 



(63) 



(64) 



again in agreement with our discussion about dominant paths. Furthermore this shows that for large /i and thus 
high filling or density, half integer values of the surface charge are thermodynamically selected, and consequently the 
applied voltage vs. the average surface charge as a function of applied voltage should exhibit plateaus at half integer 
values of Q. 

Building on these arguments, we can evaluate the value of /1 where the transition occurs. Let us compare the cases 
9 = and 9 = 0.5, choosing only the paths of maximal weight. When 9 = 0, the path of maximal weight stays at 
X n = and gives Xo(0 = 0) ~ 1. When 9 = 0.5, the path of maximal weight oscillates freely between X n = 1/2 
and X n = —1/2, leading to the eigenvalue Xo(9 = 0.5) ~ exp(— 7/8)(/i + 1). The half-integer value is preferred when 
A (6» = 0.5) > A (9 = 0), i.e. when 



fi > exp(7/8) - 1. 



(65) 



This approximate result will be compared with numerical computations later. 

In the fixed Q ensemble our approximate results and Eq. (|22p yield the following results for the voltage drop in the 
strong coupling and thermodynamic limit 



(Av) Q = M7I? for 7 6> 2 > ln(/i) 

= M-y(d - -sign(0)) for j9 2 < ln(jii). 



(66) 
(67) 



Note that when 9 is small and M is finite, the boundary terms coming from Si n a will dominate and the potential 
drop will only be weakly dependent on the system size. Interestingly we see that this extensive contribution to the 
free energy oscillates. In the limit where there are few ions (thus small fi) we see that the potential drop has the same 
sign as 9. However, when there are many ions (large /1) the potential drop has a sign opposite to that of 9. 

We emphasize that these results depend on the assumption that 7 is small, for instance at high temperature. The 
oscillatory effects seen here should weaken as the temperature is increased, however their underlying presence should 
remain visible. We also note that they should vanish when 9 is zero (i.e. for integer Q) but should be maximal when 
6 = ±1/2. 

In the constant voltage ensemble the system can minimize its bulk electrostatic energy by tuning the surface charge 
to be and integer, and thus 9 = 0; this observation thus predicts the appearance of a plateau like behavior in the 
charge as a function of the applied potential drop Av. The effective action in the constant V ensemble is then 



S*=- 7 n 2 +2(ln(^) - l)Q + QAv 

n=0,Q 



(2Q + 1)(Q + 1)Q + 2(ln(ju) - 1)Q + QAv. 



(68) 



In the limit of large Av, minimization of the above action then leads to the scaling law (Q)av 
the capacitance behaves as 



1 



CAv 



2 7 |Aw| 



I Aw 1 2 /j and thus 



(69) 



in accordance with mean field theory predictions 17] 
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V. NUMERICAL COMPUTATION 

We have, in principle, a formal solution to the thermodynamics of the LCF capacitor and explicit analytic results 
can be obtained in certain regimes. In general, however, the partition function needs to be computed numerically. 
This is quite straightforward as matrix elements K^ i decay exponentially when k and I are large. It is nevertheless 
important to keep matrix elements around n, m zero, as these are necessary in order to interact with the boundary 
term. We must thus keep modes k e [—(N — l)/2 — Q, (N + l)/2 — Q], where the N is an odd integer. The criteria 
that the matrix elements neglected by this truncation are small, means that we must choose exp(— 7fc 2 /4) to be 
small when k = ±(JV — l)/2 — Q. Practically the numerical evaluation is extremely fast and the dependence on the 
truncation value N can be eliminated by increasing N until no difference in the numerical results is seen. 

In the following numerical results, almost all the parameters are set to 1, except if the effect of their variation is 
studied. We thus have a — 1, A— lq = l, j3 = l, e = l, fi = 1. With this choice of parameters, taking /i = 1 
corresponds roughly to half- filling) . 

A. Free energy for charged boundaries 

In Fig. [5] we plot the reduced grand potential ujq as a function of the size of the system for two values of the 
charge fixed at the boundaries, for the exact solution and the mean field approximation. The dimensionless pressure 
is defined as p M+ i = ojm — um+i- We observe that 

• The free energy quickly turns linear. This asymptotic regime is reached for M ~ 2Q, implying that the structures 
of the surface layers do not significantly change once the system contains enough charge to screen the two charged 
boundary plates. In this regime of large plate separations, the pressure converges to a constant, which is the 
bulk pressure, pb = limM^oo Pm+X/2- We can then define the free enthalpy, 

g = u + Mp h , (70) 

which is presented in Fig. [3] its derivative with respect to plate separation gives the disjoining pressure, 
Pd =p-pb- 

• The bulk pressure depends on the boundary charge with periodicity one (in accordance with our analysis above). 
However the dependence is rather weak and can only be seen in systems containing more than 10 sites. 

• The interaction between the boundary plates is attractive when the separation is small, becoming repulsive 
when the size becomes larger than ~ 2Q. This is physically easy to explain: at large separations, the charges 
are completely screened by the ions. However, as in dilute electrolyte systems, there is a slight excess osmotic 
pressure due to the screening layers that are drawn into the system which gives a (weak) positive pressure at 
large plate separations. On the other hand, when the plates are close to each other, there is not enough space 
for ions to enter and screen them, they thus attract each other simply because of their opposite charge. The 
disjoining pressure goes to zero at large distances. 

• We find that that having an even number of sites is thermodynamically favorable with respect to an odd number 
of sites, this behavior is smoothed out as the size increases and is not seen in the mean field approximation. The 
physical interpretation is very simple: systems with an odd number of sites cannot be completely filled, because 
of the electroneutrality condition, that induces a high free energy cost when is large (/x = 100 in these plots). 
When the system becomes large, this effect is not important, because full configurations do not dominate the 
partition function. 

• On the mean field level, the general shape of the grand potential remains unchanged, but the oscillations between 
even and odd numbers of sites disappears. The bulk pressure is also different: it is slightly higher in the mean 
field approximation. 

We note here that the bulk pressure is a fundamental quantity that can be computed directly from the largest 
eigenvalue of the transfer matrix K, 



p h = log(Ao). 



(71) 



12 




FIG. 2: Dimensionless grand potential as a function of the system size for 7 = 1 and /j, = 100. 




FIG. 3: Dimensionless free enthalpy as a function of the system size for 7 = 1 and /1 = 100. 



B. Average potential drop for fixed surface charges 

Here we compute numerically (Au) for fixed charge Q for a large system (10 4 for the exact computation and 10 3 
for the mean field approximation). Our results are shown in Fig. 01 

The exact result can be split into two parts: a monotonic part almost equal to the mean field result, and an 
oscillating part. As predicted by our analytical treatment the oscillating part is extensive. The mean field result and 
the monotonic part of the exact result do not depend on the size of the system if it is large enough. 

The mean field result is again simple to explain: the surface charges are screened by a layer of cations and anions 
localized near the electrodes, in which the voltage drop takes place. These layers screen the surface charges so that 
the charge due to each plate seen by the bulk is zero, and thus the bulk electric field vanishes. This explains the 
absence of an extensive part in the potential drop. The oscillating, extensive part in the exact result means that the 
charges are not perfectly screened. As Q varies, the effective charge Q c g{Q), corresponding to the imperfect screening, 
oscillates with period 1, as predicted by our analytical results above. 
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FIG. 4: Average potential drop as a function of the imposed charge for 7 = 1 and fi = 1 and system size 10 4 
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FIG. 5: Left: Average surface charge as a function of the imposed voltage, for 7 = 1 and /1 = 1 with system size 10 4 . Right: 
Average surface charge as a function of the imposed voltage, for 7 = 1 and fi — 0.5 with system size 10 4 . 



C. Average charge for fixed potential difference 

The average boundary charge is shown as a function of the applied voltage in Fig. [5j The exact result increases 
in the same way as the mean field result, but when the voltage becomes large, plateaus are observed in the average 
charge, located at half-integer or integer values, (Fig. [5]), depending on the parameters 7 and fx. The plateaus show 
up only for large systems and they are smoothed out for small systems. 

The change in the plateau locations is due to the fact that integer surface charges are more stable at low fj,, whereas 
half-integer charges are favored at high /j, as we have shown in the large 7 limit in Eqs. (|62[) and (|64|) , which show 
explicitly that the largest eigenvalue of the transfer matrix, and thus the bulk pressure, depends on Q. We have 
shown the difference of the bulk pressure between Q = and Q = 0.5 in Fig. ©, and compared it to the approximate 
transition line (|65|) . The surface charge Q = is favored for small values of \i (dilute regime) and the surface charge 
Q = 0.5 is favored for large values of /1 (dense regime). This change in the plateau positions has an interesting effect 
on the capacitance near the point of zero charge and is addressed in the next section. 



14 




D. Capacitance 

We computed the capacitance as a function of the voltage for a large system (M = 10 4 ) and for different values of 
the parameters fi and 7. Two phenomena appear in the behavior of the capacitance: 

• the capacitance can show multiple oscillations which increase in number as the coupling 7 is increased and as 
the applied voltage Aw is increased. 

• These oscillations seem to be superimposed upon a trend predicted by mean-field theory with underlying bell 
and camel-like shapes, the former with a maximum at zero applied voltage and the latter exhibiting a minimum. 

Examples are shown on Fig. [7] We see that the peaks appear at low temperatures (high 7), and the camel- like 
peak appears at low densities (low fx) in accordance with mean- field theory [13] ■ The peaks are related to the plateaus 
in the average charge, so they are reduced when the system gets smaller. 

The exact results are compared to the mean field results for the bell and camel (Fig. [Sj) shaped capacitance. The 
mean-field capacitance exhibits the same behavior, without the peaks, which is coherent with the absence of plateaus. 
They also strongly resemble the capacitance appearing in jl7j ]. We also plot the capacitance curve when fi increases, 
for 7 = 1, on Fig. |H1 which shows the transition between the camel shape and the bell shape. 

We now discuss the implication of the change in the plateau positions. When they are located at integer charges, 
the PZC, Q = 0, corresponds to a plateau; the surface charge will thus respond very slowly to changes in the imposed 
voltage and the capacitance is expected to be small. On the other hand, when the plateaus are located at half integer 
charges, the PZC lies between two plateaus and the capacitance is expected to be large. In the thermodynamic limit, 
when the system size goes to infinity, we thus may even expect a sharp phase transition. We plot the capacitance 
at the PZC as a function of /1 in Fig. [TUJ we indeed see a phase transition, which becomes sharper as temperature 
decreases. 

A final point that should be mentioned is that the capacitance in the exact solution can exceed the Helmholtz 
capacitance Cjj- 
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FIG. 7: Left: Capacitance as a function of the voltage drop for fj, = 1: appearance of the peaks as 7 increases for a bell shaped 
capacitance. Right: Capacitance as a function of the voltage drop for /1 = 0.03: appearance of the peaks as 7 increases for a 
camel shaped capacitance. 
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FIG. 8: Left: Capacitance as a function of the voltage drop for fi — 1 and 7=1. Right: Capacitance as a function of the 
voltage drop for /1 = 0.03 and 7 = 0.3. 



E. Charge density 

The charge density versus the position is plotted on Fig. [TT] we show only the charge density close to the left 
boundary, the same profile with opposite charge is found close to the right boundary. It shows a counter-ion layer 
with a size corresponding to the charge Q. 

If we increase the density (// > 1) for a non integer charge Q, the charge density shows oscillations, as is seen on 
Fig. Q21 corresponding to the ion layers as observed in [22]. The phenomenology exhibited is: 

• The amplitude of the oscillations varies continuously with the charge: it is maximal for half-integer charges and 
zero for integer charges, which is physical because an integer charge can be screened perfectly, cf. Fig. [T3] 

• The mean field result does not show oscillations (as can be seen in [13, H3]) : they are thus a basic outcome of 
the exact result. 
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FIG. 10: PZC capacitance as a function of /i for different values of 7, for an infinite system. 



• We show these oscillations on a larger length scale in Fig. [T2l The oscillations are clamped with a characteristic 
length £, which increases with /i, and we notice that £ ~ Actually, the two highest eigenvalues of the transfer 
matrix K, Xq and Ai, define a correlation length via £ = (ln|Ao/Ai|) . An analytic computation for three 
Fourier modes and Q = 0, and for two modes and Q = 0.5, indeed gives £ ~ fj,. 

• The physical idea is that each boundary (left and right) creates such oscillations. On the plots described above, 
these oscillations sum up, because there is an even number of sites in the system. If there is an odd number of 
sites, the sum of the oscillations coming from the left and right boundary is destructive, as shown in Fig. [T5] 
This effect is obviously irrelevant for large (experimental) systems where M>(. 

We also note that if the char ge o n the boundary is small, the ionic densities do not saturate, which is not in 
agreement with what is found in [271 ] on the continuous mean field level. 

It can also be seen in the exact result that the surface charges are not necessarily perfectly screened (on average) by 
the ions. This can be established by integrating the average charge in the fluid over the first half of the space (which 
contains the left layer) which does not always give —Q. The average screening is perfect only when Q is integer or 
half-integer. 
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FIG. 11: Left: Charge density close to the left electrode (located at x — —1) as a function of the position for 7 = 1, y, = 1 and 
Q = 1: exact result (solid line) and mean field result (dashed line). Right: Charge density close to the left electrode (located 
at x — — 1) as a function of the position for 7 = 1, fj, = 1 and Q — 4.25: exact result (solid line) and mean field result (dashed 
line). 




The phenomenon of over-screening cab be seen in Fig. 1151 In the exact treatment, and for high densities, the charge 
of the first layer may be higher than the boundary charge, and thus over-screens it. However, no over-screening can 
be seen in the mean field approximation. We should also emphasize that over-screening is not a necessary condition 
to observe charge oscillations, e.g. oscillations exist for 7 = 1 and fi = 3 (cf. Fig. fl~2|) . a case where there is no 
over-screening. 



VI. NUMERICAL SIMULATIONS 



Because of the intricacies observed in the behavior of this system we deem it advisable to check the theoretical 
results against simulation. In this section we present a selection of simulation results which cover the different aspects 
of the ID LCF. The simulation is performed directly with the spin formulation and so the comparison serves also as 
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FIG. 13: Exact result for the mean charge density as a function of the position for 7 = 1 and /i = 100, for different values of 
the boundary charge. 




FIG. 14: Exact result for the mean charge density as a function of the position for - 
of sites (M = 81). 



1000 and Q — 0.5, for an odd number 



a check of the Hubbard-Stratonovich transformation to the formulation in terms of the field potential 4>. In all cases 
agreement between theory and experiment is absolutely excellent. 

We can consider both the canonical and grand-canonical ensembles but concentrate here on the grand-canonical 
ensemble which is appropriate to the scalar field reformulation. We simulate both for the constant Q and constant 
Av ensembles which are related by the usual Legendre transformation. From earlier, on a lattice with M sites, we 
have 

M-x 

i,j=0 

m = /3ff +7«Q+^7Q 2 (Af + l), (72) 
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FIG. 15: Mean charge density po of the first layer as a function of the surface charge Q for 7 = 1 and fi G {3, 10, 100}: exact 
result (solid line) and mean field approximation (dashed line). 



where u — X)f=o 1 For ^ e fi xe( i Q grand ensemble with chemical potential /i the partition function is 



^expU/^-logM^Tl^l 2 ) 

{S p } \ i ) 



E Q = 2_,exp( -PH-logn^Stl* 1 , (73) 



and for the fixed At; grand ensemble the partition function is 

H A „ = J dQE Q exp(QAv) 



E 

{S P } 



where Av is the fixed potential difference between the plates in the appropriate units. 

The update is done using the Metropolis algorithm. We select a pair of neighboring spins 5j, Sj+i at random. If 
\Si + 5,4-1 1 = 2 then no update is possible because of charge conservation and we select a new pair. In the case 
\S{ + 5i4.i| < 2 we suggest an update of the pair which respects charge conservation but is otherwise chosen with 
equal probability; this preserves detailed balance. We choose this probability to obtain an efficient acceptance rate. 
The suggested update is accepted/rejected using the standard Metropolis procedure. One lattice update consists of 
M pair updates. We typically perform of order 10 5 — 10 6 lattice updates to establish equilibrium and then measure 
the operator under consideration every subsequent lattice update and compute the average at the end of the run. We 
typical use 10 5 — 10 6 measurement lattice updates. All calculations were performed on multiple processor machines, 
each processor of which was used to carry out an independent simulation of the kind described here. The grand 
average over the independent simulations was taken and the error computed from the variance of the distribution of 
the individual processor averages. We used both a 12 core intel desktop computer and also an NVIDIA Tcsla GPU 
with 512 cores. Results from both these computers agree very well indeed. 

In Fig[TB] we show (Q) versus Av in the Av ensemble for 7 = 1, /i = 1, M = 128. We have 

w>-*=T"«*"W$- (75) 

The graph in the figure shows that theory and simulation agree very well. On the same figure we also show (Av) 
versus Q in the Q ensemble for 7 = 1, [i = 1, M = 1024. Here we have 

(A«)=-Aio g H = 7 ((M + 1)Q -<«)). (76) 
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FIG. 16: Left: Theory and simulation for mean charge (Q) versus At; for 7 = l,p = 1 on a lattice with M = 128. Right: 
Theory and simulation for mean voltage difference (Aw) versus Q for 7 = 1, p = 1 on a lattice with M — 1024. 



Again the agreement of simulation with theory is extremely good. 

In Fig[T7]we show the mean charge density p versus x the distance in lattice units, a, from the left hand plate which 
carries charge Q — 0.5 and is located at x = — 1. This is in the Q ensemble for 7 = 1,/z = 1,M = 1024. In this case 

P=(S X ), x = 0,1,2,... . (77) 

The layering of charge is clearly seen and the simulation faithfully reproduces the theory. Mild over-screening at x = 
can also be seen where p — —0.53 which over-screens the charge of Q = 0.5 on the plate at x — — 1. 

In Fig. [18] we show the capacitance Cav versus applied voltage Av in the fixed Av ensemble for M = 128 and the 
respective cases 7 = l,p = 1 and 7 = 0.3, p — 0.03. C& v is time-consuming to calculate since it is given in terms of 
the connected two point function (u 2 ) c : 

d 2 1 (u 2 ) 

Cav = d(A^ log Eav = t(mTT) + (MTIF (78) 

Connected two point functions are generally hard to calculate since they are a small difference between two large 
quantities for which we have separate statistical estimates. The cancelation gives an estimate that is a small number 
but with an error commensurate with that of each term separately. Our simulation results took well over 24 hours to 
complete which is why we were limited to the relatively small lattice size M = 128. The results again agree well with 
the theory curve. 

We conclude by stating that we can successfully simulate the model formulated in terms of lattice spins and that the 
results agree very well with the formulation in terms of the potential field which is given by the Hubbard-Stratonovich 
transformation. The agreement between theory and simulation also verifies the Fourier method used in the exact 
theoretical analysis. 



VII. CONCLUSION 



We have studied a ID lattice Coulomb fluid model in a capacitor configuration. The model can be solved exactly, 
allowing us to compute all the relevant quantities, notably the pressure, the capacitance and the charge distribution. 
Even though the system is one dimensional and it would thus be difficult to make quantitive comparison with experi- 
mental systems, we believe the results are pertinent to the field of ionic liquid capacitors for two main reasons. First 
it allows a comparison of an exact result with the much used mean-field theory. It thus reveals which characteristics 
of the exact solution can be captured by the mean-field theory. It seems that while the mean field theory successfully 
predicts overall trends and orders of magnitude, it fails to pick up on oscillatory behavior in charge distributions and 
some rather exotic behavior in the charge- voltage capacitor characteristics of the model. 
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FIG. 17: Theory and simulation in the fixed Q ensemble for the mean charge density p versus distance x from the left-hand 
plate located at x = —1 for y = 1, fl = 10, Q = —0.5 on a lattice with M = 1024. Note that there is mild over-screening at 
x = 0. 
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FIG. 18: Left: Theory and simulation for capacitance Cav versus Av in the fixed Av ensemble for y = 1, fi = 1 on lattice with 
M = 128. Right: Theory and simulation for capacitance Cav versus Av in the fixed Av ensemble for y = 0.3, p, — 0.03 on 
lattice with M = 128. 



The study suggests that the simple LCF model can lead to a variety of phenomena which are somehow washed out by 
the mean field approach, charge distribution oscillations, over-screening and very exotic behavior of the capacitance. 
Because the behavior predicted by our exact solution was so unexpected we have verified our analytical predictions 
by comparing them with numerical simulations of the ID model and have shown them to be accurate. Another point 
which should be made is that in contrast to [lj], H3 we developed a discrete mean-field theory on the lattice rather 
than taking a continuum limit. This is an important point as we show that even the discrete mean- field theory does 
not reproduce the various oscillatory phenomena mentioned above, which are thus not a simply consequence of a 
discrete lattice and the charge oscillations that we see are therefore due to correlations rather than being an artifact 
due to the introduction of a lattice. 

In concluding, we should also emphasize that the ID model has certain pathologies that are not present in a three 
dimensional system, notably the lack of screening of the surface charges in systems with large plate separations. If we 
would like to make a link with experimental systems it is clear that the simple picture of charged sheets would need 
to be modified. The state space of the charge on the sheets here is rather cartoon-like with just three values and just 
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two bare weights for the presence of the charges ±1 or (given by the fugacity /i). An extension of this model would 
be to introduce a larger set of dynamical values for the charges on the sheets and it may even be possible to find an 
effective one dimensional lattice model from a complete three dimensional model by exploiting the layering parallel 
to the plates that should be induced by an applied voltage, perhaps via renormalization group type arguments. 

We finally note that the results of numerical simulations agree with those of the theoretical calculations and therefore 
confirm the equivalence of the spin and scalar field formulations as well as the use of the Fourier method to carry out 
the exact theoretical calculation. 
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